{-# LANGUAGE TypeSynonymInstances, FlexibleInstances #-}
module Geometry.Intersection ( Scale (..)
                             , Rotate (..)
                             , Translate (..)
                             , Line
                             , Vector
                             , Polygon
                             , Point
                             , Angle
                             , Radius
                             , Direction (..)
                             , Circle (..)
                             , intersectionCircleWithVector
                             , intersectionCircleWithPolygon
                             , visibilityZones
                             , pointZones
                             , VisibilityZones (..)
                             , shiftVisibilityZones
                             , toPolar
                             , toAngle 
                             , inPolygon
                             , inCircle
                             , module Data.VectorSpace
                             , module Data.AffineSpace ) where

import Data.VectorSpace
import Data.AffineSpace
import Data.List (sortBy, maximumBy)
import qualified Data.IntMap as M

type Radius = Double
type Angle = Double
type Point = (Double, Double)
type Polygon = [Point]
type Vector = (Double, Double)
type Line = (Point,Point)

data Direction = In | Out deriving (Show)
data Circle = Circle Point Radius deriving (Show)

class Scale a where
  scale :: Double -> a -> a
  
instance Scale Point where
  scale s (x,y) = (x*s, y*s)

instance Scale p => Scale [p] where
  scale s = map (scale s)

instance Scale Circle where
  scale s (Circle p r) = Circle (scale s p) (r*s)
  
instance Scale (Point -> a) where
  scale s g x = g (x ^/ s)

class Translate a where
  translate :: Vector -> a -> a
  
instance Translate Point where
  translate v p = p .+^ v
  
instance Translate p => Translate [p] where
  translate v = map $ translate v
  
instance Translate Circle where
  translate v (Circle p r) = Circle (translate v p) r
  
instance Translate (Point -> a) where
  translate v g x = g (x .-^ v)

class Rotate a where
  rotate :: Double -> a -> a
  
instance Rotate Point where
  rotate a (x,y) = let x1 = x*cos a - y*sin a
                       y1 = x*sin a + y*cos a
                   in (x1,y1)

instance Rotate p => Rotate [p] where
  rotate a = map $ rotate a

instance Rotate (Point -> a) where
  rotate a g x = g $ rotate (-a) x

sides :: [a] -> [(a, a)]
sides ps = zip ps (tail ps ++ [head ps])

shift :: Angle -> Angle -> Angle
shift a b = r ((a+b)/2/pi) * 2*pi where
  r x = x - fromIntegral (floor x)
  
data QuadraticRoot = DoubleRoot Double
                   | TwoRoot (Double, Double)
                   | ComplexRoot

quadraticEquation a b c = 
  let d = b^2 - 4*a*c
      d1 = sqrt d
      x1 = (-b-d1)/2/a
      x2 = (-b+d1)/2/a
      r | d == 0 = DoubleRoot x1
        | d >= 0 = TwoRoot (x1,x2)
        | otherwise = ComplexRoot
  in r

data CircleVectorIntersection = OUT Point
                              | IN Point
                              | Tangent Point deriving (Show)

intersectionCircleWithVector :: Circle -> (Point, Vector) -> [CircleVectorIntersection]
intersectionCircleWithVector (Circle p r) (p0, v) = 
  let a = p0 .-. p
      root = quadraticEquation (v<.>v) (2*(a<.>v)) (a<.>a - r^2)
      h t = p0 .+^ v^*t
      assert t = t >= 0 && t <= 1
  in case root of
    DoubleRoot t    -> if assert t
                       then [Tangent (h t)]
                       else []
    TwoRoot (t0,t1) -> map (\(s,t) -> case s of 
                               '-' -> IN (h t)
                               '+' -> OUT (h t)) $ 
                       filter (assert . snd) [('-', t0), ('+', t1)] 
    ComplexRoot -> []
    
intersectionCircleWithPolygon :: Circle -> Polygon -> [CircleVectorIntersection]
intersectionCircleWithPolygon c = concatMap f . sides where
  f (p0,p1) = intersectionCircleWithVector c (p0, p1 .-. p0) 
                                               
check :: [(Bool, b)] -> Bool
check xs = (l `mod` 2) == 0 && (a == xs' || b == xs')
  where l = length xs
        a = take l $ cycle [True, False]
        b = take l $ cycle [False, True]
        xs' = map fst xs

data VisibilityZones = Arc [(Double, Double)]
                     | InVisible
                     | Visible
                     deriving (Show)

shiftVisibilityZones :: Double -> VisibilityZones -> VisibilityZones
shiftVisibilityZones _ Visible = Visible 
shiftVisibilityZones _ InVisible = InVisible
shiftVisibilityZones t (Arc ts) = Arc $ map (\(a,b) -> (shift a t, shift b t)) ts

visibilityZones :: Radius -> Polygon -> VisibilityZones
visibilityZones r polygon = 
  let circle = Circle (0,0) r
      ts = sortBy (\(_,a0) (_,a1) -> compare a0 a1) $ 
           map (\(s, p) -> (s, toAngle p)) $
           map (\t -> case t of
                   OUT x -> (True, x)
                   IN  x -> (False, x)) $ 
           filter (\t -> case t of
                      OUT _ -> True
                      IN _  -> True
                      Tangent _ -> False) $ intersectionCircleWithPolygon circle polygon
  in case check ts of
    False -> InVisible
    True  -> case null ts of
      True -> case inPolygon polygon (0,0) of
        False  -> InVisible
        True   -> Visible
      False -> let ts' = if (fst . head) ts == True
                         then ts else tail ts ++ [head ts]
               in Arc $ map (\[a,b] -> (a,b)) $ group1 2 $ map snd ts'

pointZones :: Point -> Polygon -> VisibilityZones
pointZones p polygon = shiftVisibilityZones (negate (toAngle p)) $ visibilityZones (magnitude p) polygon

group1 n [] = []
group1 n xs = let (a,b) = splitAt n xs
              in a : group1 n b

inCircle :: Circle -> Point -> Bool
inCircle (Circle p r) p1 = magnitude (p1.-.p) <= r

inPolygon :: Polygon -> Point -> Bool
inPolygon ps p = and $ map (f p) $ sides ps where
  f :: Point -> (Point,Point) -> Bool
  f p (p1,p2) = 
    let a = p2 .-. p1
        b = p .-. p1
    in a `det` b >= 0

det :: Vector -> Vector -> Double
det (a,b) (c,d) = a*d - b*c

toAngle :: Point -> Angle
toAngle = snd . toPolar

toPolar :: Point -> (Radius, Angle)
toPolar (x,y) = 
  let m = magnitude (x,y)
      a = acos (x/m)
  in if y >= 0
     then (m, a)
     else (m, 2*pi - a)

rectangle :: Polygon
rectangle = scale 0.5 [(1,1), (-1,1), (-1,-1), (1,-1)]

